!> date: 2022-10-23
!> 保存为二进制 VTU 文件
module sph_save_vtk

    use seakeeping_string, only: to_string
    use seakeeping_filesystem, only: operator(.join.)
    use vtk_fortran, only: vtk_file
    use sph_region_type, only: region_t
    implicit none

contains

    !> 保存为 VTK/VTU 文件
    subroutine save_vtk(vtk, ik, region)
        use sph_io_env, only: idx, ones
        character(*), intent(in) :: vtk         !! vtk 所在文件夹
        integer, intent(in) :: ik               !! 第 ik 个时间步
        type(region_t), intent(in) :: region    !! 计算域
        integer :: stat
        type(vtk_file) :: vtu

        stat = vtu%initialize(format='binary', &
                              filename=vtk.join.'pif'//to_string(ik, 'i0')//'.vtu', &
                              mesh_topology='UnstructuredGrid')

        stat = vtu%xml_writer%write_piece(np=region%ntotal, nc=region%ntotal)
        stat = vtu%xml_writer%write_geo(np=region%ntotal, &
                                        nc=region%ntotal, &
                                        x=region%loc(:, 1), &
                                        y=region%loc(:, 2), &
                                        z=region%loc(:, 3))
        stat = vtu%xml_writer%write_connectivity(nc=region%ntotal, &
                                                 connectivity=idx(0:region%ntotal - 1), &
                                                 offset=idx(1:region%ntotal), &
                                                 cell_type=ones(:))
        stat = vtu%xml_writer%write_dataarray(location='node', action='open')
        stat = vtu%xml_writer%write_dataarray(data_name='Velocity', &
                                              x=region%vel(:, 1), &
                                              y=region%vel(:, 2), &
                                              z=region%vel(:, 3))
        stat = vtu%xml_writer%write_dataarray(data_name='Mass', &
                                              x=region%mass)
        stat = vtu%xml_writer%write_dataarray(data_name='Pressure', &
                                              x=region%p)
        stat = vtu%xml_writer%write_dataarray(data_name='Density', &
                                              x=region%rho)
        stat = vtu%xml_writer%write_dataarray(data_name='Itype', &
                                              x=region%itype)

        stat = vtu%xml_writer%write_dataarray(location='node', action='close')
        stat = vtu%xml_writer%write_piece()
        stat = vtu%finalize()

    end subroutine save_vtk

end module sph_save_vtk
